Reproduction time statistics and segregation patterns in growing populations 
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Pattern formation in microbial colonies of competing strains under purely space- limited population 
growth has recently attracted considerable research interest. We show that the reproduction time 
statistics of individuals has a significant impact on the sectoring patterns. Generalizing the standard 
Eden growth model, we introduce a simple one-parameter family of reproduction time distributions 
indexed by the variation coefficient S £ [0, 1], which includes deterministic (5 — 0) and memory- 
less exponential distribution (S = 1) as extreme cases. We present convincing numerical evidence 
and heuristic arguments that the generalized model is still in the KPZ universality class, and the 
changes in patterns are due to changing prefactors in the scaling relations, which we are able to 
predict quantitatively. At the example of Saccharomyces cerevisiae, we show that our approach 
using the variation coefficient also works for more realistic reproduction time distributions. 

PACS numbers: 87.18.Hf, 89.75.Da, 05.40.-a, 61.43.Hv 



I. INTRODUCTION 

Spatial competition is a common phenomenon in 
growth processes and can lead to interesting collective 
phenomena such as fractal geometries and pattern for- 
mation [IHII . This is relevant in biological contexts such 
as range expansions of biological species d, Q or growth 
of cells or microorganisms, as well as in social contexts 
such as the dynamics of human settlements or urbaniza- 
tion Q. These phenomena often exhibit universal fea- 
tures which do not depend on the details of the particu- 
lar application, and have been studied extensively in the 
physics literature [1&0-G3- 

Our main motivating example will be growth of mi- 
crobes in two dimensional geometries, for which recently 
there have been several quantitative studies. In general, 
the growth patterns in this area are influenced by many 
factors, such as size, shape and motility of the individual 
organism (Tlj . as well as environmental conditions such 
as distribution of resources and geometric constraints 
[HI fl3j , which in turn influence the proliferation rate or 
motility of organisms [l4|. We will focus on cases where 
active motion of the individuals can be neglected on the 
timescale of growth, which leads to static patterns and 
is also a relevant regime for range expansions. We fur- 
ther assume that there is no shortage of resources, and 
growth and competition of species is purely space lim- 
ited and spatially homogeneous. This situation can be 
studied for colonies of immotile microbial species grown 
under precisely controlled conditions on petri dish with 
hard agar and rich growth medium. 

Under these conditions one expects the colony to form 
compact Eden- type clusters [l2j , which has recently been 
shown for various species including Saccharomyces cere- 
visiae, Escherichia coli, Bacillus subtilis and Serratia 
marcescens [Til . [l5j . 

The Eden model [l6| has been introduced as a basic 
model for the growth of cell colonies. It has later been 



shown to be in the KPZ universality class 0, 0, Eli fl8l |. 
which describes the scaling properties of a large generic 
class of growth models. In recent detailed studies of 
E. coli and S. cerevisiae [pfl . ll5L [l9l I20I ] quantitative ev- 
idence for the KPZ scaling of growth patterns has been 
identified. The models used in these studies ignored all 
microscopic details reproduction, such as anisotropy of 
cells [2l[, and could therefore not explain or predict dif- 
ferences observed for different species. Nevertheless, they 
provided a good reproduction of the basic features such as 
KPZ exponents, which is a clear indication that segrega- 
tion itself is an emergent phenomenon. Fig. [1] shows dif- 
ferences in growth patterns in a circular geometry taken 
from fTo| for immotile E. coli and S. cerevisiae. For both 
species the microbial populations are made of two strains, 
which arc identical except having different fluorescent la- 
beling. Reproduction is asexual, and the fluorescent la- 
bel carries over to the offspring. At the beginning of 
the experiments the strains are well mixed, but during 
growth rough sector shaped segregated regions develop. 
The qualitative emergence of these segregation patterns 
and connections to annihilating diffusions has been stud- 
ied in [l5|, [l9|, HO, [22j , ignoring all details specific for a 
particular species. 

For S. cerevisiae the domain boundaries are less rough 
when compared to E. coli, leading to a finer pattern con- 
sisting of a larger number of sectors. In general, this 
is a consequence of differences in the mode of reproduc- 
tion and shapes of the microbes, which introduce local 
correlations that are not present in simplified models. 
In this paper we focus on the effect of time correlations 
introduced by reproduction times that are not exponen- 
tially distributed (as would be in continuous time Marko- 
vian simulations), but have a unimodal distribution with 
smaller variation coefficient. This is very relevant in most 
biological applications (see e.g. |23l - [2a |). and even in 
spatially isotropic systems the resulting temporal corre- 
lations lead to more regular growth and therefore smaller 
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FIG. 1: (Color online) Fluorescent images of colonies of (a) 
E. coli and (b) S. cerevisiae. The scaling properties of both 
patterns are believed to be in the KPZ universality class, and 
the differences are due to microscopic details of the mode of 
reproduction and shape of the micro-organisms. The images 
have been taken with permission from [la ], copyrighted (2007) 
by the National Academy of Sciences, U.S.A. 



fluctuations of the boundaries, with an effect on the pat- 
terns as seen in Fig. [2 

To systematically study these temporal correlations, 
we introduce a generic one-parameter family of reproduc- 
tion times, explained in detail in Section ITT] We establish 
that the growth clusters and competition interfaces still 
show the characteristic scaling within the KPZ univer- 
sality class, and the effect of the variation coefficient is 
present only in prefactors. We predict these effects quan- 
titatively and find good agreement with simulation data; 
these results are presented in Section Hill More realistic 
reproduction time distributions with a higher number of 
parameters are considered in Section [IV! where we show 
that to a good approximation the effects can be sum- 
marized in the variation coefficient and mapped quan- 
titatively onto our generic one-parameter family of re- 
production times. Therefore, our results are expected to 
hold quite generally for unimodal reproduction time dis- 
tributions, and the variation coefficient alone determines 
the leading order statistics of competition patterns. 



II. THE MODEL 

For regular reproduction times with small variation co- 
efficient the use of a regular lattice would lead to strong 
lattice effects that affect the shape of the growing clus- 
ter. To avoid these we use a more realistic Eden growth 
model in a continuous domain in R 2 with individuals 
modelled as circular hard-core particles with diameter 
1, since we want to study purely the effect of time cor- 
relations. This leads to generalized Eden clusters which 
are compact with an interface that is rough due to the 
stochastic growth dynamics. 

Let B(t) denote the general index set of particles p 
at time t, (x p ,y p ) <E R 2 is the position of the centre 
of particle p, and s p £ {1,2} is its type. We write 
B{t) = B\{t) U B2{t) as the union of the sets of particles 
of type 1 and 2. We also associate with each particle 
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FIG. 2: (Color online) A smaller variation coefficient S in 
reproduction times (see @ and (JSJ) leads to more regu- 
lar growth, smoother domain boundaries and finer sectors. 
Shown are simulated circular populations with (a) 5 = 1 and 
(b) 5 = 0.1. Both colonies have an initial radius of ro = 50, 
and they are grown up to simulation time t = 50 leading to 
final radii of approximately 120 (a) and 95 (b). The different 
colors denote cell types 1 and 2. 

the time it tries to reproduce next, T p > 0. Initially, T p 
are i.i.d. random variables with cumulative distribution 
F$ with parameter 5 e [0, 1], which is explained in de- 
tail below. After each reproduction T p is incremented by 
a new waiting time drawn from the same distribution. 
Note that we focus entirely on the neutral case, i.e. the 
reproduction time is independent of the type and both 
types have the same fitness. We describe the dynamics 
below in a recursive way. 

Following a successful reproduction event of particle p 
at time t = T p , a new particle with index q — \B(T p — )| + 1 
is added to the set B Sq with the same cell type s q = s p , 
such that 

B 3p (T p +) = B Sp (T p -)U{q} . (1) 

Here B(T p —) and B(T p +) denote the index set just be- 
fore and just after the reproduction event, and |-B(t)| 
denotes the size of the set B(t). The position of the new 
particle is given by 

(x q ,y q ) = (x P ,y P ) + (cose/), sine/)) , (2) 

where (f> G [0, 27r) is drawn uniformly at random. This 
is subject to a hard-core exclusion condition for circular 
particles, i.e. the Euclidean distance to all other particle 
centres has to be at least 1 , as well as to other constraints 
depending on the simulated geometry as explained be- 
low. Note that in our model the daughter cell touches 
its mother, which is often realistic but in fact not essen- 
tial, and the distance could also vary stochastically over 
a small range. The new reproduction times of mother 
and daughter are set as 

T p H- T° ld + T , T q = T° ld + T' , (3) 

where T, T' are i.i.d. reproduction time intervals with 
distribution F$. There can be variations on this where 
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mother and daughter have different reproduction times, 
which arc discussed in Section IIVI The next reproduc- 
tion event wiil then be attempted at t = min \T q : q G 
B(T p +)}. Reproduction attempts can be unsuccessfui, if 
there is no available space for the offspring due to block- 
age by other particles. In this case the attempt is aban- 
doned and T p is set to oo, as due to the immotile nature 
of the cells this particle will never be able to reproduce. 

The initial conditions for spatial coordinates and types 
depend on the situation that is modelled. In this paper 
we mostly focus on an upward growth in a strip of length 
L with periodic boundary conditions on the sides, where 
we take -B(O) = {!,... ,L} with (x p , y p ) = (p, 0), for all 
p E B(0). The initial distribution of types can be either 
regular or random depending on whether we study single 
or interacting boundaries, and will be specified later. 

In Section Hm for our main results we use reproduction 
times T distributed as 



1 - S + Exp(l/<5) , 5 G (0, 1] 



(4) 



i.e. T has an exponential distribution with a time lag 
1 — S € [0, 1) and a mean fixed to (T) = 1 for all S. 
The corresponding cumulative distribution function F$ 
is given by 



Fs(t) 





1 _ e -(t-l+5)/5 



t<l-S 
t>l-6 



(5) 



The variation coefficient of this distribution is given by 
the standard deviation divided by the mean, which turns 
out to be just 



V {T 2 ) — (T) 2 _ S _ 
(T) 1 



(6) 



With this family we can therefore study reproduction 
which is more regular then exponential with a fixed av- 
erage growth rate of unity (equivalent of setting the unit 
of time). 

For 6=1 this is a standard Eden cluster, but S < 1 
introduces time correlations. While the correlations af- 
fect the fluctuations, we present convincing evidence that 
they decay fast enough not to change the fluctuation ex- 
ponents, so the system remains in the KPZ universality 
class. Furthermore we make quantitative predictions on 
the (5-dcpcndcncc of non-universal parameters and com- 
pare them to simulations. The more synchronized growth 
leads to effects similar to the ones seen in experiments 
(Fig. [1]). To give a visual impression of the patterns 
produced by the model, we show in Fig. [5] two growth 
patterns with 5=1 and 0.1. The initial condition is a 
circle, and the types are distributed uniformly at random. 
The patterns are qualitatively similar to the experimen- 
tal ones in Fig. [TJ and more regular growth leads to a 
finer sector structure. The same effect is shown on Fig. [3] 
for the simulations in a linear geometry with periodic 
boundary conditions, which is analyzed quantitatively in 
the next Section. Smaller values of 6 also lead to more 
compact growth and smaller height values reached in the 
same time. 




FIG. 3: (Color online) Populations in a linear geometry with 
periodic boundary conditions in lateral direction with (a) S = 
1 and (b) 5 = 0.1. Both populations have lateral width L = 
300, and the colonies are grown to a simulation time t ~ 50, 
leading to heights of approximately 70 (a) and 40 (b). The 
different colors denote cell types 1 and 2. 



III. MAIN RESULTS 
A. Quantitative analysis of the colony surface 

In this Section we provide a detailed quantitative anal- 
ysis of the 5 family of models in linear geometry with 
periodic boundary conditions (see Fig. [3]) , starting with 
the dynamical scaling properties of the growth interface. 

We regularize the surface to be able to define it as a 
function of the lateral coordinate x and time t as 

y(x,t) := max{y p : p G B(t), \x p - x\ < l} . (7) 

In case of overhangs (which are very rare) we take the 
largest possible height, and due to the discrete nature of 
our model this leads to a piecewise constant function. 

The surface of a standard Eden growth cluster is known 
to be in the KPZ universality class [HI EH, i.e. a suitable 
scaling limit of y{x, t) with vanishing particle diameter 
fulfills the KPZ equation 

dty(x,t) = v + uAy(x,t) + ^{Vy{x,t)f + s/D n {x,t). 

(8) 

Here vq, of the order of unity, corresponds to the growth 
rate of the initial flat surface (related to the mean repro- 
duction rate and some geometrical effects), the surface 
tension term with v > represents surface relaxation, 
and the nonlinear term represents the lowest order con- 
tribution to lateral growth [l8|]. The fluctuations due 
to stochastic growth are described by space-time white 
noise r](x,t), which is a mean Gaussian process with 
correlations 

{ n (x,t)r,{x',t l ))=8{x-x , )8{t-t l ). (9) 

We denote the average surface height by 



(10) 



which is a monotone increasing function in t. It is also 
asymptotically linear and therefore we will later also use 
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FIG. 4: (Color online) Family- Vicsek scaling (|12[) of the sur- 
face roughness S(t). The data collapse under rescaling with 
a = 1/2 and z = 3/2 occurs in a scaling window which is 
narrower for small 5 due to intrinsic correlations. The dif- 
ferent symbols correspond to different values of 5, and the 
color represents system size, L — 1500 (green/light grey) and 
L = 4000 (blue/dark grey). The dashed lines indicate the 
expected slope ft = 1/3. The data for L = 1500 has been 
averaged over 100 samples and for L = 4000 over 30 samples. 
The error bars are comparable to the size of the symbols. 



h as a proxy for time. The <5-depcndence of the average 
growth velocity of height as seen in Fig. [3] does not lead 
to leading order contributions to the statistical properties 
of the surface or the structure of sectoring patterns. 

The roughness of the surface is given by the root mean 
squared displacement of the surface height as a function 
of t @, 03, defined as 



s{t) = (l j~ (vM-Ht)) 2 dx) 1/2 . (11) 

The main properties of the surface y(x, t) can then be 
characterized by the Family- Vicsek scaling relation of the 
roughness 

S(t) = L a f(t/L z ) , (12) 
where the scaling function f(u) has the property 

J 3 u < 1 



/(«) 



1 



u > 1 



(13) 



Such a scaling behaviour has been shown for many dis- 
crete models inclu ding ballistic deposition and continuum 
growth 0, EH EH, l26l [27J , and holds also for other uni- 
versality classes such as Edward Wilkinson. For the KPZ 
class in 1 + 1 dimensions the saturated interface roughness 
exponent is a = 1/2, the growth exponent is ft = 1/3, 
and the dynamic exponent is z = a/ (3 = 3/2. 

Fig. 0] shows a data collapse for the roughness S(t) for 
two system sizes, and for a number of different values 
of 5. As 5 gets smaller, the surface becomes less rough 



FIG. 5: (Color online) The height-height correlation function 
C(l,t) for L = 4000 at t = 11000 for various values of 5. The 
data has been averaged over 30 samples, and the error bars 
are comparable to the size of the symbols. The dashed lines 
indicate the expected slope 1/2. 



due to a more synchronized growth. The dashed lines 
indicate the power law growth with exponent /3 = 1/3 in 
the scaling window. This window ends around t/L z « 1 
due to finite size effects, where the lateral correlation 
length reaches the system size and the surface fluctua- 
tions saturate. For small t the system exhibits a transient 
behaviour before entering the KPZ scaling due to local 
correlations resulting from the non-zero particle size and 
stochastic growth rules. As we quantify later, these cor- 
relations are much higher for more synchronized growth 
at small S, which leads to a significant increase in the 
transient regime. The transient time scale is indepen- 
dent of system size and vanishes in the scaling limit, so 
that the length of the KPZ scaling window increases with 
L. This behaviour can be observed in Fig. U where for 
the smallest value (5 = 0.2 the scaling regime is still hard 
to identify for the accessible system sizes. 

Another characteristic quantit y is the height-height 
correlation function defined as pj.l28l [29| 

C(l, t) = (j J\y{x, t) - y(x + I, t)) 2 dx) 1/2 . (14) 

For a KPZ surface in 1 + 1 dimensions this obeys the 
scaling behaviour 



C(M) ~i(l) 2/3 (A*)V3 



(15) 



where £n (t) is defined to be the lateral correlation length 
scale and takes the form 0, [2i| |3(| 



£|i(t) - (L»/2^) 1 / 3 (At) 2 / 3 



(16) 



A detailed computation can be found in Appendix [Bj For 
small values C(l,t) grows as a power- law with I, and when 
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FIG. 6: (Color online) Dependence of the KPZ parameters 
D/(2v) on 8. Data are obtained from ()15|1 by fitting the pref- 
actor of the power law in Fig.0 using that the proportionality 
constant is very close to 1 (cf. derivation (|B11[) in the Ap- 
pendix). The data are in good agreement with the prediction 
pT) with fitted parameters e « 0.42 and D/(2v)(5 = 1) » 1.1. 



Z exceeds the lateral correlation length £y(i) it reaches a 
value that depends on the time t and the parameters 
of ©. This is shown in Fig. [SJ where C(l,t) is plotted 
for various values of 5, and the data agree well with the 
exponent a — 1/2 for the KPZ class indicated by dashed 
lines. 

The time correlations introduced by the partial syn- 
chronization can be estimated by considering a chain of 
N growth events where each particle is the direct de- 
scendant of the previous one. Each added particle corre- 
sponds to a height change Aj/j, and has an associated 
waiting time T{ with distribution During time t 

there are N(t) growth events, and since the average re- 
production time is 1 with variance 5 2 , we have (N(t)) ~ t 
and var(JV(t)) ~ 5 2 t. The height of the last particle is 

VN(t) = £f W A Vh leading to 
var(j/ JV(t) ) = (Ai M ) 2 var(JV(t)) + (JV(t)) var(A yi ) . (17) 

The terms in this expression correspond to two sources 
of uncertainty: (i) due to the randomness in Ti the num- 
ber of growth events vary with var(iV(i)), and (ii) the 
individual height increments are random with var(Ay.;). 
This leads to 

var( 2/7V(t) )«t(Ay l ) 2 (5 2 + e 2 ) , (18) 



where e = \/ var( Aj/j) / (Ajfj) denotes the variation coeffi- 
cient of the height fluctuations due to geometric effects. 

We define the correlation time r as the amount of time 
by which the uncertainty of the height of the chain be- 
comes comparable to one particle diameter, var(yjv(r)) = 
0(1). Since (Ayj) is largely independent of 6 (cf. Ap- 
pendix |X| , the time correlation induces a fixed intrinsic 



vertical correlation length 



in the model. This correlation length reduces fluctuations 
and leads to an increase in the saturation time t sat of 
the system, namely t sat /r ~ L z , a modification of the 
usual relation with the system size L. Analogous to the 
standard derivation of the time-dependence of the lateral 
correlation length [|| , this leads to 

£||(i)~(t/r) 1/z . (20) 

Together with (|16p from the behaviour of the correlation 
length, we expect 

D/(2v)~(5 2 + e 2 ) 2 , (21) 

since A turns out to be largely independent of 5. This 
is shown to be in very good agreement with the data 
in Fig. El for fitted values of e and a prefactor. The fit 
value for e and the ratio D/(2v) for S = 1 (the usual 
Eden model) are compatible with simple theoretical ar- 
guments (see Appendix [A"]l . So the very basic argument 
above to derive an intrinsic vertical correlation length 
explains the (5-dependence of the surface properties very 
well. Measuring height in this intrinsic length scale, we 
observe a standard KPZ behaviour with critical expo- 
nents being unchanged, since the intrinsic correlations 
are short range (i.e. decay exponentially on the scale r) . 
This is in contrast to effects of long-range correlations 
where the exponents typically change, see e. g. studies 
with long-range temporally correlated noise [3ll - l33l ] or 
memory and delay effects using fractional time deriva- 
tives and integral/delay equations j33 - l36j . 

B. Domain boundaries 

In this section we derive the superdiffusive behaviour 
of the domain boundaries between the species from the 
scaling properties of the interface. Since the boundaries 
grow locally perpendicular to the rough surface, they are 
expected to be superdiffusive, which has been shown for 
a solid on solid growth model in (37} and has been ob- 
served in [Tol l for experimental data. In order to con- 
firm this quantitatively for our model, we perform simu- 
lations with initial conditions Bi(0) = {1, . . . , [L/2]} and 
-82(0) = {[L/2] + 1, . . . , L}, i.e. the initial types are all 
red on the left half and all green on the right half of 
the linear system. Therefore we have two sector bound- 
aries X\ and X2 with initial positions A^i(O) = 1/2 and 
X 2 (0) = [L/2] + 1/2. After growing the whole cluster, 
we define the boundary as a function of height via the 
leftmost particle in a strip of width 2 and medium height 
h: 

X x {h) = min {x p + 1/2 : \y p - h\ < 1, p G B±} 
X 2 (h) =mm{x p + l/2:\y p -h\<l,peB 2 } ,(22) 
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FIG. 7: (Color online) Scaling behaviour of the mean square 
displacement M(h) (|24jl . The system size is L = 1000 and 
the data is averaged over 500 samples and the error bars are 
comparable to the size of the symbols, (a) Data collapse of 
the normalized quantity M(h)/ag as a function of height h 
for several values of S. The values in the normalization erf are 
taken from the best fit shown as full line in (b). Each curve 
follows a power law with exponent 4/3, the line corresponding 
to /i 4//3 is shown as comparison, (b) The prefactor erf, where 
the data are best fits according to (|24p . The solid line used 
for the collapse in (a) follows the prediction (5 2 + e 2 ) 4 / 3 with 
fitted e = 0.40, which is compatible with the fit in Fig. [6] 



where we take the periodic boundary conditions into ac- 
count. The simulations are performed on a system of size 
L = 1000, and run until a time of t = 2000, which is well 
before the expected time of annihilation, which is of or- 
der X 3 / 2 proportional to the saturation timcscale in the 
KPZ class. Therefore we can treat the sector boundaries 
as two independent realizations of the boundary process 
(X(h) :h>0). 

As has been noted already in [33], this process is ex- 
pected to follow the same scaling as the lateral correlation 



length. For the mean square displacement 

v 2 



M(h) :={(X(h)-X(0)y 



(23) 



we therefore get with ([TfJ]) and (|2"Tj) . using the linear re- 
lationship between h and i, 



M(h) Kajh 



2 u2H 



em- 



(24) 



Here er 2 cx 



(S 2 + e 2 ) 4 ' 3 and the Hurst exponent is 
H = 2/3, which quantifies the supcrdiffusivc scaling of 
the mean square displacement (|23[) . This prediction is in 
very good agreement with data for the scaling of M(h) 
and its prefactor as presented in Fig. [7J and the fit value 
for e 2 is consistent with the one in Fig. [51 As before, 
for D I (2z/) the independence is absorbed by the prefac- 
tor, and the power law exponent 4/3 for M(h) remains 
unchanged from standard KPZ behaviour studied also in 

We can further investigate the law of the process 
(X(h) : h > 0). The data presented in Fig. E^a) clearly 
support that X{h) is a Gaussian process. A fractional 
Brownian motion with stationary increments seems to 
be a natural model for the X(h) in the KPZ scaling 
window. This is confirmed by the behaviour of the cor- 
relation function (X(h + Ah)X(h)), which is shown in 
Fig. EJb) for various 6 and two values of the lag Ah > 0. 
For a fractional Brownian motion with mean square dis- 
placement ([Ml) we expect 



(X(h+Ah)X(h)) a -j-{{h+Ah) 2H - 



■2H 



-\Ah\ 



2H\ 



(25) 

for all Ah > and h > sufficiently large to have no 
effects from the flat initial condition. For simplicity we 
have assumed here that X(0) = 0. 

This is in good agreement with the data, and we con- 
clude that the sector boundaries can be modeled by frac- 
tional Brownian motions with supcrdiffusivc Hurst expo- 
nent H = 2/3 and a <5-dcpcndcnt prefactor as (JM]). We 
note that the exponent H = 2/3 has also been observed 
in experiments [15j. 



C. Sector patterns 

In [l!| , and also [2fj| . [22| under the assumption of diffu- 
sive scaling, it was shown how the understanding of the 
single boundary dynamics leads to a prediction for sector 
statistics for well-mixed initial conditions. In this section 
we shortly review this approach and show that it carries 
over straight away to systems with 5 < 1. The sector 
boundaries Xi{h) interact by diffusion limited annihila- 
tion which drives a coarsening process, as can be seen in 
Fig. [3] for two linear populations with different values of 
5. Both systems have the same initial condition with a 
flat line of particles of independently chosen types, and 
the finer coarsening patterns for smaller values of 5 re- 
sult from the reduced boundary fluctuations due to the 
prefactor as ([2~4"ll . 
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FIG. 8: (Color online) The sector boundary X(h) behaves like 
a fractional Brownian motion, (a) The standardized prob- 
ability density function (pdf) of X(h) as a function of the 
rescaled argument x / (osh 2 ^) for different heights h and val- 
ues of 5. The black solid parabola is the pdf of a standard 
Gaussian on a logarithmic scale, (b) The covariance function 
(X(h + Ah)X(h)\ shows the behaviour (|25p . which is plotted 
as the solid black curve. After rescaling we get a data collapse 
as a function of h/Ah, which agrees well with the prediction 
if h is sufficiently large and the flat boundary conditions be- 
come irrelevant. Data are averages over 1000 realizations and 
the error bars are comparable to the size of the symbols. 



FIG. 9: (Color online) The average number of sector bound- 
aries (N(h)) follow a power law (|26|l with exponent —2/3, 
which is indicated by the full line. The data are plotted for a 
system size L — 1500 and various values of 5 (see legend), and 
collapse on the function h~ 2 ^ s when rescaled by L/ ^k<j 2 . 
Data are averages over 500 realizations and the error bars are 
comparable to the size of the symbols. 



system size L in the rescaling so that rescaled quantities 
are of order 1, and all data collapse on the function /i~ 2 / 3 
without prefactor. 

Using (|26[) . we can predict the expected number of 
sector boundaries at the final height in the simulations 
shown in Figure [3] For 5=1, the final height is h as 70 
leading to (N(h)) w 7.6, and for 5 = 0.1, h « 40 with 
(N(h)) ps 32. These numbers are compatible with the 
simulation samples shown which have 6 and 34 sector 
boundaries remaining, respectively. 

In general, diffusion limited annihilation is very well 
understood, and there are exact formulas also for higher 
order correlation functions [4(j, which can be derived 
from the behaviour of a single boundary (j2"4")) . This 
demonstrates that the behaviour of populations is funda- 
mentally the same for all values of d and characterized by 
the KPZ universality class, and the observed difference 
in coarsening patterns can be explained by the functional 
behaviour of the prefactors. 



Let N(h) be the number of sector boundaries at height 
h > as defined in (l22l) . For systems of diffusion limited 
annihilation (38l . l39j it is known that N(h) is inversely 
proportional to the root mean square displacement, and 
decays according to 



(N(h)) 



1 



1 



^4ttM(/i) crs 



(26) 



This prediction is confirmed in Fig. [3J where we plot 
(N(h)) for various S, and obtain a data collapse by mul- 
tiplying the data with ^ A-ko^/L, [39[. We include the 



IV. REALISTIC REPRODUCTION TIMES 

In this section we study the effect of more realistic re- 
production time distributions on the sectoring patterns, 
and how they can be effectively described by the previ- 
ous 5-dcpcndcnt family of distributions in terms of their 
variation coefficient. As an example, we focus on S. cere- 
visiae, which is one of the species included in Hal and 
for which reproduction time statistics is available [23l - [25j 
by the use of time lapsed microscopy. S. cerevisiae cells 
have largely isotropic shape so that spatial correlations 
during growth should be minimal, fitting the assump- 
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FIG. 10: (Color online) Comparison of realistic reproduction 
times with the S model, (a) The probability density functions 
of reproduction times of mother cells T r (full orange line/light 
grey) and daughter cells T m + T r (dashed orange line/light 
grey) with normalized mean compared to T from Q with 
corresponding S (blue/dark grey), (b) The prefactor of the 
mean square displacement <rf as introduced in (|24|l and Fig. [Jj 
The data correspond to reproduction times T r for all cells 
(denoted M), T m + T T for all cells (denoted D) and the most 
realistic mixed model (denoted M and D) as explained in the 
text. All these cases are consistent with previous results from 
Fig. El 



tions of our previous model. However, the results of this 
section cannot be applied directly to quantitatively pre- 
dict the patterns in Fig. [I] since the variation coefficients 
under the experimental conditions in [l5[ are not known 
to us. 

When yeast cells divide, the mother cell forms a bud on 
its surface which separates from the mother after growth 
to become a daughter cell. The mother can then im- 
mediately restart this reproduction process, whereas the 
daughter cell has to grow to a certain size in order to 
be classed as a mother and to be able to reproduce. We 
denote this time to maturity by T m and the reproduction 
time of (mother) cells by T r . 

The results in [23T - I251 . [ill ] suggest that Gamma distri- 
butions arc a reasonable fit for the statistics for Tm and 



T r is distributed as po + Gamma(/3i, p-i) 



(27) 



with delay po > 0. The parameters pi,p2 denote the 
shape and scale of the Gamma distribution, which has a 
probability density function 



f M _, P i-l eX P(-V^) 



t > o 



The time to maturity is 

T m distributed as Gamma(/?3, p^), 



(28) 



and in [25j data are presented for which the parameters 
can be fitted to po ~ 1-0, p\ « 1.7, p2 ~ 0.51, P3 ~ 9 and 
Pi w 0.3. The unit of p$, pi and p^ are hours and pi, p% 
are dimcnsionlcss numbers. 

The random variables T m and T r may be assumed to be 
independent and the time until a newly born daughter 
cell can reproduce is distributed as the sum T m +T r . Note 
that the expected value of reproduction times (T r ) = 
P0+P1P2 = 1-86 is smaller than that for times to maturity 
(T m ) = /J3/J4 = 2.52 but of the same order. The real 
time scale for these numbers is hours, but we are only 
interested in the shape of the distributions rescaled to 
mean 1 like our previous model. 

The distribution (g]) of 5-depcndent reproduction times 
can be written as T = 1 — S + Gamma(l, 5), since expo- 
nentials are a particular case of Gamma random variables 
with shape parameter 1. The reproduction time T r of 
mother and T m + T r of daughter cells are also unimodal 
with a delay, and very similar in shape to T in our model. 
This can be seen in Fig. HUT a). where we plot the prob- 
ability densities rcnormalized to mean 1. Analogous to 
([6]), we can compute the variation coefficients of T r and 
(T m + T r ), which turn out to be 0.356 and 0.244, respec- 
tively. To confirm that the behaviour of sector bound- 
aries can be well predicted by the variation coefficient, we 
present data of three simulations in Fig. fTOl b): one with 
reproduction times T r for mother and T<i + T r for daugh- 
ter cells as explained above, one with T r for all cells, and 
one with T r + T m for all cells. The mean square displace- 
ment M(h) for these models also shows a power law with 
exponent 4/3 analogous to Fig. [71 and the prefactors as 
match well with our simplified model. 

To estimate the variation coefficient in the model with 
mother and daughter cells, we measure the fraction of 
reproduction events of daughter cells to be pd = 0.88, 
and p m = 0.12 for mother cells. The reproduction time 
of the union of mother and daughter cells is then taken 
as 

T distributed as 9(T m + T r ) + (1 - Q)T r , (29) 

where the independent Bernoulli variable = Be(pd) G 
{0, 1} indicates reproduction of a daughter. The varia- 
tion coefficient of T turns out to be 0.322. In all three 
combinations of realistic reproduction times we find that 
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the generic family of Fs introduced in (U) provides a good 
approximation for the properties of domain boundaries in 
simulations. We expect this method of mapping realistic 
reproduction time distributions to our ^-dependent fam- 
ily to hold for a large class of microbial species which 
have similar unimodal distributions. 



V. CONCLUSION 

We have introduced a generalization of the Eden 
growth model with competing species, regarding the re- 
production time statistics of the individuals. This is 
highly relevant in biological growth phenomena, and can 
have significant influence on the sectoring patterns ob- 
served e.g. in microbial colonies. Although growth of im- 
motilc microbial species is the prime example, our results 
also apply to more general phenomena of space limited 
growth with inheritance, where the entities have a com- 
plex internal structure that leads to non-exponential re- 
production times, such as colonization/range expansions 
or epidemic spreading of different virus strands. Our 
main result is that, as long as the reproduction time 
statistics have finite variation coefficients, the induced 
correlations are local and the macroscopic behaviour of 
the system is well described by the KPZ universality 
class. The dependence of the relevant parameters in that 
macroscopic description on the variation coefficient (a 
microscopic property of the system) is well understood 
by simple heuristic arguments, which we support with 
detailed numerical evidence. 

Figures [2] and [3] illustrate that changes in the varia- 
tion coefficient S of reproduction times lead to significant 
changes in the competition growth patterns in our model, 
and we are able to quantitatively predict this dependence. 
We studied the effects of reproduction time statistics in 
a generalized Eden model, isolated from other influences 
such as shape of the organisms or correlations between 
mother and daughter cells which might be relevant in real 
applications. In that sense our results are of a theoretical 
nature. However, they indicate that the variation coeffi- 
cient of reproduction times can have a strong influence on 
observed competition patterns. This coefficient has been 
measured for various species in the literature, where it 
is found that it depends on experimental conditions such 
as type of strain, concentration of nutrients, temperature 
etc. [42h44( . For example, it was found that for S. cere- 
visiae the coefficient for mother cells can vary between 
S re 0.12 - 0.38 and for daughter cells 5 re 0.19 - 0.28 
depending on concentrations of guanidinc hydrochloride. 
It has also been observed that S can be as small as 0.047 
for these yeast type organisms [45j . For E. coli values 
of S re 0.32 - 0.51 have been observed in 0, HI El, 
which is larger, and compatible with the observations 
in Fig. [TJ But for the experimental conditions in [l5[ 
with pattern growth the coefficient has not been deter- 
mined and therefore the results in this paper cannot be 
readily applied to explain the differences in competition 



patterns between S. cerevisiae and E. coli. In partic- 
ular, the latter have anisotropic rod shape which has 
probably a strong influence and is currently under in- 
vestigation in the group of O. Hallatschek. Another rod 
shaped bacterium, Pseudomonas aeruginosa, has varia- 
tion coefficient 6 re 0.14- 0.2 [4lLli5|. This bacterium 
along with E. coli belongs to the gram-negative bacteria 
family. Despite obvious similarities between P. aerugi- 
nosa and E. coli in the shape of their cells, their colonics 
display morphological differences (48|, which fit qualita- 
tively into our results, and would be another interesting 
example for a quantitative study. 

In general, it is an interesting question if the sim- 
ple mechanism of time correlations due to reproduc- 
tion time statistics with variable variation coefficients is 
sufficient to quantitatively explain sectoring patterns in 
real experiments. We are currently investigating this for 
S. cerevisiae which are approximately of isotropic shape, 
in order to quantify the influence of other factors on 
growth patterns, such as correlations between mother 
and daughter cells. For example, an effective attraction 
between cells which is often observed in the growth of 
microbial colonies would influence the growth directions, 
and further smoothen the surface and the fluctuations of 
sector boundaries. For future research, it should also be 
possible to describe spatial effects due to non-isotropic 
particle shapes with the methods used in this paper. 
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Appendix A: Effect of geometric disorder 

The squared variation coefficient e 2 in p^|) due to ge- 
ometric disorder has been consistently fitted to values 
around 0.4 with our data in Section lTHl This value is com- 
patible with the following very simple argument. Con- 
sider a single growth event around an isolated spherical 
particle with diameter 1, with direction a chosen uni- 
formly in a cone with opening angle ir/2 around the ver- 
tical axis. This leads to (Ay.;) = Jl^ 2 C0Sa ^ ~ 0-64 
and 

e 2 re ( / cos 2 a— - (Ay,) 2 ) /(Ay,) 2 re 0.23 , (Al) 

which is of the same order as the fitted values. Choosing 
only a slightly larger opening angle 0.557T of the cone 
leads to e 2 re 0.39 and (Ay t ) re 0.57. These are in good 
agreement with the fitted values and with measurements 
of (Ayi) (not shown). The latter show some dependence 
on S, related also to the compactness of growth as seen 
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in Figs. [3] and [3J but this does not contribute to our 
results on a significant level so we ignore this dependence. 
Actual growth events in the simulation are of course often 
obstructed by neighbouring particles, but the right order 
of magnitude of the parameters can be explained by the 
basic argument above. 



Appendix B: Deriving the correlation function C(l,t) 

We use the mode coupling method [2|| , in order to find 
an exact analytical expression of the correlation function 
Eq. ([T4"]) as shown in Eq. ([15]) . The idea of the mode 
coupling approximation is that properties of solutions of 
the KPZ equation ([8]) may be derived by first considering 
the linear Edwards- Wilkinson equation J4|| for A = 0. 
We further consider the co-moving frame, so that vq = 0, 
and the equation then reads 



d t y(x, t) = vAy(x, t) + VDr)(x, t). (Bl) 
We denote by 



y(k,t) 



dx y(x, t)e 



-ikx 



the Fourier transform of the function y(x,t). The evolu- 
tion of the function y(k,t) satisfies 



d t y(k, t) = -vk 2 y(k, t) + VD fj(k, t). (B2) 

Here f](k, t) is the spatial Fourier transform of the white 
noise rj{x, t), where ?)(fc, t) has a mean with correlations 



(fj(k, t)f)(k', £')) = S{k + k')5{t - t'). (B3) 



A formal solution of (|B2[) can be obtained, and after in- 
verse Fourier transform this leads to 



the correlation function ([Til) of the Edwards- Wilkinson 
equation 



D 



C(l,tf = — / dfcfc- 2 (l-cos(fc/))[l-e 



-2uk 2 U 



(B7) 



In order to compute the correlation function for the 
KPZ equation ([5]) we substitute length scale dependent 
parameters D(k) and v(k) into (|B7I) . which are obtained 
from the renormalization group flow equations 0, 
In d = 2 dimensions these are given by 

v{k) = Vl [{l-a B ) + a B /k] 1 ' 2 , 
D(k) = D 1 [(l-a B ) + a B /k} 1 / 2 , (B8) 



and A(fc) = Ai, where 



a B 



X\Dx 



Here (\i,v\,Di) are the parameters for k = 1 where no 
renormalization has taken place. Plugging this into (|B7|) 
and only considering the most dominant terms, we obtain 



T) /*°° 

C (l,t) 2 = ^ / dfcfc-^l-cosOfcOMl-e- 5 ^ 2 '] 



where B = 2iy 1 a 1 J 2 = ^±\ y fD 1 /2v 1 . 
If we take t — > oo in Eq. (]B9p we get 



(B9) 



C(l,t) 2 ^ Dl 



dk k- 2 (l- cos(W)) = — I . (BIO) 

With (|B8[) D/v — Di/vi is independent of the scale k, 
and thus 

C(Z,t)« (§;l) 1/2 for i<C||(*)- (BH) 



/OO pt 
dke lkx / dsf ] {k,s)e- vk2{t ^ s) . (B4) 
-oo JO 

The correlation function C{l,t) defined in ([14]) can be 
represented as 

C(l,t) 2 = 2[ dx(y(x : t) 2 -y(l+x,t)y(x,t)) . (B5) 
Jo 

Using the solution (|B4[) we can compute 

'dx{y(l+x,t)y(x,t)) = £-p k C -^[l-e- 2 » k2t ] , 

(B6) 

where we have used that the Fourier transform is even 
in k. Taken together, this leads to an expression for 



For finite time, numerical integration of (]B9I) in the large 
I limit gives 



lim C(l,t) 2 ps 2.7— (Bt) 



2/3 



Together with ([Blip and the definition ([14]) of the corre- 
lation length this leads to 



lim C(l,t) ps U.A x 2 1 ^(—Y /3 7r~ 5 / 3 (Xt) 2 ^) 
z-s-oo V \2v ) J 



and 



^(t)«5.4x2 1 /3g) 1/ V 5 / 3 (A0 



5/3/^^/3 



1/2 

(B12) 
(B13) 
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